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A method based on a set of new theorems for the analysis of multichannel time series is 
described, based on precise Fourier transform and coherence analysis of the restored 
signals from a detailed set of frequency components. Magnetic field recordings of 
spontaneous and evoked activity by means of magnetic encephalography demonstrated 
that multichannel precise Fourier spectrum contains a very large set of harmonics with 
high coherence. The inverse problem can be solved with great precision based on 
coherent harmonics, so the technique is a promising platform of general analysis in brain 
imaging. The analysis method makes it possible to reconstruct sites and timing of electrical 
activity generated by both spontaneous and evoked brain function at different depths in 
the brain in the millisecond time range. 

Keywords: multichannel signal analysis, frequency coherence, extraction of partial spectra, time series 
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INTRODUCTION 

Modern scientific studies are performed by means of new pow- 
erful equipment, generating large amounts of detailed data. 
Magnetic encephalography (MEG) provides an example of 
a foremost biological technology, comparable with the most 
sophisticated physical devices. Magnetic encephalographs reg- 
ister magnetic field in hundreds of channels with sampling 
frequency up to several thousand Hertz. Typical 5 min exper- 
iment on the 275 channel device with sampling rate 1200 Hz, 
provides 100 million field values, so the problem of data anal- 
ysis appears a pressing challenge in the MEG technique. Many 
approaches are used to solve various scientific and diagnostic 
problems of encephalography. Fourier analysis in many imple- 
mentations can be called the oldest and the most popular of 
methods used for the brain data analysis (Dietsch, 1932; Jansen 
et al., 1981; Halliday et al., 1995). Through the whole history 
of this method it was connected with difficulties of calcula- 
tions, so the development of the Fast Fourier Transform (FFT) 
(Cooley and Tukey, 1965) dramatically advanced the application 
of the Fourier analysis in many fields, including brain research 
(Miyashita et al., 2003). Regularization of the FFT was performed 
in multitaper method (Thomson, 1982; Percival and Walden, 
1993), implemented in the studies of the evoked and sponta- 
neous activity (Llinas et al., 1999; Mitra and Pesaran, 1999). 
In the quantitative electric and magnetic encephalography stud- 
ies, trying to study patterns of the Fourier coefficients, rough 
spectral analysis is combined with statistical analysis of coher- 
ence between channels or independent components (Tauscher 
et al, 1998; Weiss and Rappelsberger, 2000; Jarvis and Mitra, 
2001; Makeig et al, 2002; Delorme and Makeig, 2004; Garcia- 
Rill et al, 2008; Muthuraman et al, 2010). Usually in appli- 
cations of the Fourier analysis to brain studies the spectra are 



calculated in short (<10s) time windows, based on the well- 
known property of instability of the brain processes (Llinas, 
2001). 

Typically brain studies register activity in many channels 
simultaneously for protracted time periods (up to tens of min- 
utes in hundreds of channels). Those registered data are usually 
processed with two important methodological weaknesses: First 
weakness is that in time dependence analysis the methods are 
applied, which were developed for the solitary time series, mul- 
tichannel recordings are implemented mainly to attempt inverse 
problem solutions. The second weakness lies in the usage of short 
time windows (less than 10 s), decreasing the resolution of the 
Fourier transform. 

These procedural limitations relate to the fact that only 
descriptive conclusions can be drown from such studies. Thus, 
it is often stated, that the particular pathology or evoked activ- 
ity results in the spectral changes of some kind or another. Such 
qualitative approach diminishes a detailed analysis quite drasti- 
cally, while it can be useful from a general diagnostic perspective 
of view or to study some general cognitive features. At the same 
time that approach to data analysis loses practically all experimen- 
tal information, and narrows encephalography to general sets of 
particular observations. 

Recently the method of precise frequency-pattern analysis to 
decompose complex systems into functionally invariant entities 
was proposed (Llinas and Ustinin, 2012). The method makes it 
possible to address general spectra to the partial spectra of stable 
functional entities and to restore their time series. The method is 
based on the complete utilization of the long time series, while 
the multichannel nature of the data is also completely taken into 
account, making it possible to implement detailed reconstruction 
of neuronal circuit activity. 
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METHODS 

In multichannel recording of the brain activity, such as those from 
MEG, a magnetic field event is recorded by many channels at 
discrete time moments, providing sets of experimental vectors 
{bk}, where fc is number of channel, the Z-th component of vector 
bj; = bk(l) is the result of field measurement at the time moment 
tj, 1= 1, ...,L. 

The first step in the method proposed is the interpolation of 
the experimental data in every channel, providing the continuous 
function 

hit) = F(b h t), (1) 

where F is a function, corresponding to the particular method of 
interpolation (Boyd, 2001). The linear and spline interpolation 
methods are used in our study with satisfactory results. 

After interpolation the multichannel recorded signal is 
described by the set of functions {Bk(t)}, where t e [0, T] , T = 
ti — t\, T is the whole time of measurement, k = 1 , . . . , K, K is 
maximal number of channel. 

The multichannel precise Fourier transform calculates a set of 
spectra for interpolated functions {Bk(t)}: 



2 f T ~ 2 f T - 

«o* = ~ / B k (t)dt, a nk = - B k (t)cos(2nv n t)dt, 
1 Jo 1 Jo 



Kk=7z\ B k (t)sm{2nv n t)dt, 
1 Jo 



(2) 



where a ok , a nk , b nk are Fourier coefficients for the frequency v n 
in the channel number k, and n = 1, . . . , N, N = v max T, where 
Umax is the highest desirable frequency. 

The term "Precise" is used in three different senses here and is 
achieved by three distinct steps: 

1. Precise calculation of the Fourier integrals. Gaussian quadra- 
ture formulas are used to calculate integrals on any interval [0, 
T] , in the registration scale. 

2. Building all spectra for the total registration time T, as opposed 
to methods using moving or fractional window. The step in 
frequency is equal to Av = v„ — v„_ i = ^, thus frequency 
resolution is determined by the recording time. 

3. Tuning of the frequency grid by cutting the interval of integra- 
tion T to build the optimal approximation of the frequency 
selected. Tuning can be performed by little changes of the 
integration time T. 

The method can be implemented without interpolation, if 
Fourier integrals can be calculated with a required accuracy on 
the experimental set of time points f/, I = 1, . . . ,L and if tuning 
of the frequency grid is not necessary. 

This precise transform leads to an accurate and reversible 
representation of time data in the frequency domain for each 
channel. As for the space domain, "space" is determined by the 
simultaneous registration by multiple channels having differ- 
ent positions with respect to the source. That is, if an accurate 
representation of time series for all channels is used, spatial 
characteristics of the signal can also be determined accurately. 



Given a precise multichannel spectra it is possible to perform 
the inverse Fourier transform using: 



B k (t) = — + a nk cos(2jtv n t) + b nk sin(27rv„f), 



v„ = — , N = v max T 



where aok, a nk , b nk are Fourier coefficients, found in (2). 
This formula can also be written as 



(3) 



"Ok 



Bh ' ' ' = if + ^ P " k Sin( - 2jZVnt + tpnk ">' V " = ~f ' N = V max T ' ( 4 ) 
n=l 



where p nk = ^ ' a\ k + b 2 nk , cp„ k = atan2(a„jfc, b nk ). 

The transform (3) or (4) allows the possibility of implement- 
ing precise filtering, including or eliminating any selected set of 
frequencies when restoring the signal. 

We propose to study the detailed frequency structure of the 
brain, restoring multichannel signal at every frequency and ana- 
lyzing the patterns obtained. 

The multichannel signal is restored at particular frequency in 
all channels: 



B„k(t) = p„k sin(27r V„f + (p nk ), 



(5) 



where t e [0, T Vn j, k= 1, . . . , K and T Vn = -j- is the period of 
this frequency. 

The summary instantaneous power produced by all channels 
at the frequency v n will be: 



Pn(t)=J2 B lk(t)- 



The proximity of phases cp n k in different channels can be char- 
acterized by the value of empirical one-frequency coherence: 



Ci f = 1 



min (p„(t)) 
t6[o,r„„] 

max (p„(t))' 
te o,r v „ 



(6) 



where min and max are calculated at the period T Vn . Possible 
values of the coherence lay between 0 and 1: C\f € (0, 1]. 
The physical sense of C\f follows from formula (5). If all 
channels have equal phases <p nk = (p„ at the frequency v„, 
then Cif is equal to 1. If phases in different channels are 
distributed uniformly and amplitudes are equal, then Cif 
approaches to zero when maximal number of channels -K" is 
growing. 

The pattern of magnetic field at the time moment t is deter- 
mined by relation between values of the induction in different 
channels and by their average energy. Relative values make it pos- 
sible to determine the spatial structure of the source from the 
inverse problem solution, and this structure is the same for the 
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same relative values of the channels. If <p„k = <p n , then formula 
(5) can be written as 

Bnk(t) = Pnk sin(2jr v„t + (p n ) = P nk p„ sin(2jr v n t + <p n ), (7) 

where p n = JzLiPtk and Pnk = P f n ■ 

From the formula (7) it follows that relative values of chan- 
nels P n k are independent of time. The inverse problem solutions, 

determined by the normalized pattern P„k, have the same spatial 
structures for any moment of the restored time. The amplitude of 
the source is determined by p„ sin (2nv n t + ip n ) — common for all 
channels, meaning that this source is oscillating as a whole at the 
frequency v n . 

The following theorems have been proved. 

Coherence Theorem 1. The equality of phases in all channels 
is a necessary and sufficient condition for normalized pattern 
invariance through reconstructed time. 

Conclusion 1. If for particular frequency phases are equal in 
all channels, then the spatial structure of the source at this fre- 
quency can be found by the solution of inverse problem for the 

pattern p n k. 

Coherence Theorem 2. The equality of phases in all chan- 
nels is necessary and sufficient condition for the equal- 
ity of the empirical one-frequency coherence to 1, C\f = 1. 
This theorem provides a directly calculable feature to esti- 
mate the proximity of phases in all channels at any 
frequency v n . 

Coherence Theorem 3. The time course of the magnetic field 
source, having arbitrary spatial structure, can be restored from 
the partial Fourier spectrum. This partial spectrum consists of 
the frequencies with C\t equal or close to 1, having the same 
normalized pattern. Spatial structure of the source can be found 
from this pattern. 

Consider the equivalent current dipole (ECD), characterized with 
two vectors: ro — is the location of the dipole and Q — is a dipo- 
lar moment. The model of ECD in spherical conductor (Sarvas, 
1987) is used to calculate the magnetic induction registered by 
sensor, having the location r and direction n: 

B(ro, Q) = -^t (CF(Q x r 0 ) - (Q x r 0 , r) VF) , n) , (8) 

where F = a(ar+r 2 - (r 0 , r)), Vi 7 = (a 2 r _1 + cT 1 (a, r) + 2a 
+ 2r) r - (a + 2r + cT 1 (a, r)) r 0 , 

a = r — ro, a = \a\ , r = |r| , |n| = 1, (Iq = 4tv ■ 10~ 7 . 

It can be shown, that magnetic induction depends linearly on the 
dipole moment and can be written as: 

B (r 0 , Q) = ((r 0 x n)F - (VP, n) (r 0 x r)) , q) = (K,Q) , 

(9) 

where K = (F (r 0 x n) - (VF, n) (r 0 x r)). 



From the formula (9) and the principle of superposition it fol- 
lows, that the induction, measured by the sensor number k from 
/ dipoles can be written as (Hamalainen et al., 1993): 

; 

B k = (% Q;) • (10) 

;'=i 

Consider a coherent system, consisting of J dipolar sources, hav- 
ing similar time dependencies: Qj = Cj Q; Q(t), where Q(t) is a 
function of time, common for all dipoles, Cj gives the force of 

dipole number j, Q, is a unitary vector, giving the direction of 
this dipole. 

The formula (10) now can be written as 

; 

ha) = Q(t)J2 p k,cj, (id 

where Pkj = (^-kj, QjJ ■ The lead field matrix element Pkj is given 
by the sensing character of the sensor number k in relation to the 
source number j. After the summation 

; 

Pk = Y, P kj c )> d 2 ) 

hit) = Q(t)p k . d3) 

After the precise Fourier transform it follows from the for- 
mula (13) that for every frequency the restored signal in the k-th 
channel will be 

B„k(t) = PkPn sm(27tv n t + <p n ), (14) 

where p n sin (2nv n t + ip n ) = Q„(f) is the n-th Fourier compo- 
nent of the function Q{t). 

Formula (14) can be written in the form (7): 

B„k(t) = PPkPn sin(27rv„f + <p n ), (15) 

where P = J^L iH and h = §- 

The instantaneous power will be: p n {t) = P 1 p\ sin 2 (27rv„f + 
q> n ), and it has minimum = 0 in the period, so C\f is equal to 1 for 
every frequency of the function Q(f). From (15) it also follows, 
that normalized patterns at all frequencies of the function Q(f) 
will be identical. This makes possible to extract the partial spec- 
trum, corresponding to this source, from the full multichannel 
spectrum. This is implemented by selecting frequencies with high 
coherence, having equal or similar normalized patterns. From 
this partial spectrum the time course Q(f) can be restored by the 
inverse Fourier transform. 

These theoretical considerations are the foundation for the 
reconstruction of time courses of static functional entities (neural 
circuits, or sources), performing detailed frequency analysis and 
studying the similarity of the patterns with high coherence. Also 
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spatial structure of the sources at separate frequencies with high 
coherence can be restored, leading to the total decomposition of 
the brain activity. 

The algorithm of mass precise frequency-pattern analysis can 
be summarized as: 

(1) Precise Fourier Transform of the multichannel signal. 

(2) Inverse Fourier Transform — restoration of the signal at each 
frequency. 

(3) If the coherence at the particular frequency is close to 1, 
then use the pattern and frequency as elementary coherent 
oscillation. 

(4) If the restored signal consists of several phase-shifted coher- 
ent oscillations, then extract those oscillations. 

After the fourth step the initial multichannel signal will be 
represented as a sum of elementary coherent oscillations. 
Each elementary oscillation has distinct frequency, constant 
pattern and being produced by the functional entity having 
constant spatial structure. The set of elementary coherent 
oscillations is unique for the subject and for the particular 
experiment. 

(5) Split the spectrum to the assembly of partial spectra, based 
on the extraction of frequencies with similar normalized 
patterns. 

(6) Solve inverse problem for normalized patterns in order 
to determine spatial structure of the functionally invariant 
entities. 

(7) Restore the time courses from partial spectra of functionally 
invariant entities. 

After the seventh step all brain activity recorded is decomposed 
to the set of activities of functional entities, with known spatial 
structures and time courses. 

The method can be called "Frequency-Pattern Functional 
Tomography," because it reveals the structure and function of the 
brain under study in the particular experiment. 

EXPERIMENTAL RESULTS 

The method proposed makes it possible to perform detailed study 
of the brain structure and function by means of multichannel 
measurements, such as magnetic or electric encephalography. 
Note, that Theorems 1-3 make sense only if the coherence of 
multichannel oscillations is high at the frequency under study. 
In order to estimate the applicability of the method to the real 
data, MEG experimental data sets for 19 persons (control sub- 
jects and patients) were processed. Nine data sets were obtained 
with a 148-channel magnetometer Magnes 2500 WH, and 10 data 
sets were obtained with a 275-channel gradiometer. Both sys- 
tems were installed in the Bellevue Hospital in the Center for 
Neuromagnetism of New York University School of Medicine. 
These experiments were performed during research studies of 
human spontaneous and evoked activity, including control sub- 
jects and patients with various disorders (Llinas et al, 1999; 
Ustinin et al., 2010). The NYU Institutional Review Board and 
Bellevue Hospital Research Protocol Review Group approved the 
study and an informed written consent was obtained from all 
subjects before the MEG recording. 



As an example of the Precise Fourier analysis, consider the 
auditory experiment JG03_01. The MEG data were obtained with 
a 275-channel synthetic third order gradiometer (VSM MedTech 
LTD) at Bellevue Hospital in the Center for Neuromagnetism at 
the Department of Neuroscience and Physiology of the New York 
University School of Medicine. The auditory stimulus was a 2 ms 
click, presented at 14 Hz into the left ear of a healthy subject, and 
the recording was implemented. 

The precise multichannel Fourier spectrum was calculated, 
using the whole registration time (~300s) and sampling rate 
1200 Hz. Figure 1 illustrates the multichannel spectrum in the 
frequency range 1.5-50 Hz. This spectrum was calculated with 
Av = v n — v„_i = 3i3qHz and contains 15,000 frequencies in 
275 channels, plotted simultaneously. It can be concluded, that 
this plot gives only general impression about the spectrum, and 
detailed quantitative analysis is necessary. Such precise analysis is 
illustrated in Figures 2-4. 

The tuning of the frequency grid was done by minor cutting 
the time of registration T to build an optimal approximation of 
the basic stimulus frequency and its harmonics (Figure 2). This 
tuning increases the signal to noise ratio by the order of magni- 
tude in power, giving the possibility to solve the inverse problem 
with high precision. 

A large scale in frequency was utilized which produces a pre- 
cise spectrum obtained that is fractioned into separate frequencies 
with different patterns of Fourier amplitudes. 

The fragment of the spectrum in narrow frequency band 
shown in Figure 3, is close to the second harmonic of the stim- 
ulus. The normalized patterns P„k (7) of the restored magnetic 
fields are also presented for five neighboring frequencies. Also 
illustrated in Figure 3 is the extraction of the response to stim- 
ulus from spontaneous activity. The frequency peak at the second 
harmonic of the stimulus produces normalized pattern, corre- 
sponding to auditory cortical activity. The frequency to analyze 
was selected from the precise Fourier transform of the stimu- 
lus time course. Following the inverse Fourier transform of this 
frequency, highly coherent oscillation presents second harmonic 
of the brain response to auditory stimulus. The structure of the 
functional entity, generating this response, can be found from the 
inverse problem solution (Figure 4). 

In order to study general properties of the restored multichan- 
nel Fourier components, the statistical analysis of the coherence 
distribution was performed (Figure 5). For every frequency, the 
inverse Fourier transform (5) was performed and the empirical 
coherence (6) was calculated. Also the summary energy for every 
frequency was calculated. 

It was found that, in this experiment, mean coherence was 
0.51. Moreover 4.42% of frequencies have coherence greater than 
0.8 and 0.25% of frequencies have coherence greater than 0.9. 
(Figures 5E,F). This illustrates the fact, that the highest amplitude 
frequencies have the greatest coherence. 

For example, 4.42% of all frequencies, with coherence greater 
than 0.8, represent 11.26% of all energy, while 0.25% of fre- 
quencies with coherence greater than 0.9 produce 1.87% of total 
energy. Statistical properties of the restored multichannel one- 
frequency oscillations are illustrated by Figure 5. The conclusion 
is that there are many frequencies with high coherence, providing 
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FIGURE 2 | Tuning of the precise multichannel Fourier spectrum for the 
stimulus frequency in auditory MEG experiment JG03 01. 

(A) Multichannel spectrum close to the second harmonic of the stimulus, 
calculated using T = 300 s. (B) The same spectrum, calculated using 



T = 300 — 0.017 s. It can be concluded, that a negligible lost of information 
significantly increases the signal to noise ratio (S/N), leading to much higher 
coherence. More precise solution of the inverse problem can be obtained 
after such optimization. 



the direct information about the functional entities, generating 
those frequencies. 

Precise multichannel Fourier spectra were calculated for all 
19 data sets, using (2). For all frequencies of the every spectrum 
the multichannel inverse Fourier transform (5) was performed, 
the empirical coherence C\f (6) was calculated, and the analysis 
described in Figures 1-5 was performed. The results were simi- 
lar to those described in Figures 1-5, with the average coherence 
0.523. 

It was found, that: 

(1) Precise multichannel spectra of the MEG data parcel to 
oscillations with different patterns. 



(2) Both in spontaneous and evoked experimental spec- 
tra there are many frequencies with coherence, 
close to 1. 

(3) MEG sensory evoked activity allow attaining data with high 
coherence at the selected stimulus frequencies, which dis- 
close the brain functional structure generating particular 
responses. 

The conclusion is that Theorems 1-3 can be applied to 
the experimental data in multichannel magnetic and electric 
encephalography, giving a direct nonparametric method to reveal 
functional entities in the human brain, oscillating as whole 
systems. 
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In some cases, when applying the above method to real data, 
the value of empirical coherence is less than 1 . 
This may occur for various reasons, namely: 

1. Non-correlated noise, produced by the system under study, 
including sensors noise. This noise can be reduced through 
experimental design, for example, increasing the time of mea- 
surement T. 



35- 



□ □□□□ 




27.99 



FIGURE 3 | Detail of the Figure 2B, close to the maximum of spectrum 
at 28 Hz, (the 2nd harmonic of the stimulus). Five normalized patterns 
P„k (7) of the restored MEG, corresponding to neighboring frequencies, are 
shown. 



2. Activity of different non-correlated sources, falling at the same 
frequency band v n ± ^y. 

This is typical for most methods utilizing Fourier analysis, espe- 
cially when moving or fractional windows are utilized. The precise 
Fourier transform can address the issue imaging activity from 
different non-correlated sources by either increasing the record- 
ing time T or/and by tuning of the frequency grid to an exact 
frequency. 

3. Activity of several coherent sources, shifted in phase, having 
exactly the same frequency v n and also falling at the same 



frequency band v n ± 



2T- 



In the third case, simultaneous activity of coherent sources with 
different spatial structures can indicate functional connectivity. 
To divide different coherent processes from the restored mul- 
tichannel signal independent component analysis of patterns 
activity can be used, or those processes can be separated directly. 
This procedure leads to the extraction of patterns, produced by 
several different sources with high coherence at the same fre- 
quency. It can reveal the physiological dependence of the activity 
sources at this frequency. 

DISCUSSION 

A very large set of multichannel MEG recordings have been pro- 
cessed using precise frequency-pattern analysis. Following this 
procedure it has been found, that for many frequencies, empir- 
ical coherence (formula 6) is close to 1. As follows from the 




FIGURE 4 | The inverse problem solution reveals the spatial structure of 
the functional entity, producing the pattern of auditory response 
(Figure 3, middle pattern). This structure is well described by the two-dipole 



model. (A-D) Show the sections of the tomogram, going through the right 
dipole. (E) Shows the map of experimental magnetic field. (F) Shows the 
map of the field, produced by the inverse problem solution. 
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FIGURE 5 | Statistical properties of the restored oscillations. C lf — one 

frequency empirical coherence (6). (A) Shows relative number of 
frequencies, falling in the 5%-bin of coherence. (B) Shows share of 
frequencies, having coherence greater than chosen value. (C) Shows 
relative summary energy of frequencies, falling in the 5%-bin of coherence. 
(D) Shows share of common energy, produced by frequencies, having 
coherence greater than chosen value. (E) Shows distribution of the relative 
one-frequency energy, falling in the 5%-bin of coherence. For example, 
average energy, produced by one frequency with coherence >0.95, is 11 
times greater than average energy, produced by one frequency. (F) Shows 
folded cumulated distribution functions, calculated from (B) blue and 
(D) red. 



Coherence Theorem 1, normalized pattern for the particular fre- 
quency with high coherence is close to invariant. This means that 
the functional entity, producing this pattern, is oscillating as a sin- 
gle entity at this frequency. Thus, the activity of the elementary 
subsystem of the brain can be extracted and its spatial structure 
determined from the inverse problem solution. In the particu- 
lar case of magneto-encephalography this spatial structure is very 
robust as no arbitrary parameters were used to extract the pat- 
tern. Also in case of the MEG the inverse problem solution does 
not need subject parameters, except the head shape (Hamalainen 
etal., 1993). 

Given the above the technique described is capable of address- 
ing the structure of the brain as a set of coherent functional 
entities. 

Note, that this method can also be applied in electroen- 
cephalography. In the case of EEG the pattern will be also 
extracted correctly, but the solution of the inverse problem will 
involve complicated brain parameters, such as the spatial distri- 
bution of conductivities. Those parameters are usually unknown, 



adding the uncertainty to the inverse problem solution (Fuchs 
et al., 2001; Hallez et al., 2005; Plis et al, 2007; Barnes et al, 2008). 

Finally, our results indicate that precise spectra are sui generis 
every subject, in particular concerning spontaneous activity. We 
propose this new analysis paradigm for brain research, based on 
the calculation of precise spectra and on their storage for future 
reference concerning the development of pathological conditions, 
among other uses. The general number of functional entities in 
the particular experiment can be estimated as 5-10 thousand. 
This number in the order of magnitude is close to the num- 
ber of categories, introduced to describe cognitive processes in 
Huth et al. (2012). Functional entities, revealed by the method 
proposed (Llinas and Ustinin, 2012), correspond to emergent 
functions of neural circuits (Alivisatos et al., 2012) and the study 
of these entities can be the important component of the starting 
Brain Activity Map project. 
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